#ifndef RSIM_CHOLESKY_H
#define RSIM_CHOLESKY_H

/* -------------------------------------------------------------------------- *
* File: Cholesky.h                                                           *
* Authors: Kishor Bhalerao                                                   *
* Email : kishor8dm@gmail.com                                                *
* Contributors:                                                              *
* Email:                                                                     *
* copyright (c) 2010 Authors.                                                *
*                                                                            *
* This program is free software: you can redistribute it and/or modify it    *
* under the terms of the GNU General Public License as published by the Free *
* Software Foundation, either version 3 of the License, or any later version.*
*                                                                            *
* This program is distributed in the hope that it will be useful, but        *
* WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY *
* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License   *
* for more details. (http://www.gnu.org/licenses/)                           *
* -------------------------------------------------------------------------- */

#include "Matrix.h"
#include "LinAlgDecl.h"

#include <math.h>

namespace RSIM
{

/** 
<P>
Note: This file was obtained from Template Numerical Toolkit (TNT).
Details of TNT can be found on the website http://math.nist.gov/tnt/
This file is minimally modified to make it compatible with RSIM.
The originial file name (available for download on TNT website) 
is jama_cholesky.h.

<P>
kdb:This class does not have as much safety checks as i would like to have.
Hope to include more checks in the future versions.

<P>
For a symmetric, positive definite matrix A, this function
computes the Cholesky factorization, i.e. it computes a lower 
triangular matrix L such that A = L*L'.
If the matrix is not symmetric or positive definite, the function
computes only a partial decomposition.  This can be tested with
the is_spd() flag.

<p>Typical usage looks like:
<pre>
        Matrix A(n,n);
        Vector L;

        ... 

        Cholesky chol(A);

        if (chol.is_spd())
                L = chol.getL();
                
        else
                cout << "factorization was not complete.\n";
        </pre>
<p>
        (Adapted from JAMA, a Java Matrix Library, developed by jointly 
        by the Mathworks and NIST; see  http://math.nist.gov/javanumerics/jama).

*/

class Cholesky
{
private:
        // lower triangular factor
        Matrix L_;
        
        // 1 if matrix to be factored was SPD
        int isspd;

public:
        /**
        Constructs a lower triangular matrix L, such that L*L'= A.
        If A is not symmetric positive-definite (SPD), only a
        partial factorization is performed.  If is_spd()
        evalutate true (1) then the factorizaiton was successful.
        */
        Cholesky(const Matrix& A);
        
        /**
        @return the lower triangular factor, L, such that L*L'=A.
        */
        const Matrix& getL() const{return L_;}
        
        /**
        Solve a linear system A*x = b, using the previously computed
        cholesky factorization of A: L*L'.
        
        @param  B   A Matrix with as many rows as A and any number of columns.
        @return     x so that L*L'*x = b.  If b is nonconformat, or if A
        was not symmetric posidtive definite, a null (0x0)
        array is returned.
        */
        Vector solve(const Vector& B);
        
        
        /**	
        Solve a linear system A*X = B, using the previously computed
        cholesky factorization of A: L*L'.
        
        @param  B   A Matrix with as many rows as A and any number of columns.
        @return     X so that L*L'*X = B.  If B is nonconformat, or if A
        was not symmetric posidtive definite, a null (0x0)
        array is returned.
        */
        void solve(const Matrix& B, Matrix& A);
        
        /** 
        Solves the linear system Ax = B. CAUTION. No dimension checks anywhere in
        this function.
        */
        void solve(const Matrix& B, double *x);
        
        
        void inverse(double *x);
        
        /**
        @return 1, if original matrix to be factored was symmetric 
        positive-definite (SPD).
        */
        int is_spd() const {return isspd;}

};
} // namespace RSIM

#endif
